;+
;PURPOSE
;  convert SLUG binned histograms to non-binned stars
;SYNTAX
;	hess2stars, modname[, path=path, abslumcut=abslumcut]
;INPUTS
;	modname: the name of the SLUG model (e.g. the default
;		for SLUG is "SLUG_DEF"
;	path: set the path to the directory which has the
;		histo_modname.fits file [default to '$SLUG_DIR/output']
;	abslumcut: since this file can get quite large,
;		this allows an absolute luminosity cut on the 
;		stars which are output (in the band
;		defined by Y_AXIS_NAME on the histo_modname.fits file)
;		[defaults to 3]. Setting this to a
;		high value will dramatically increase output file size
;OUTPUT
;	the output file will be written to
;          path+'cmd_'+modname+'.fits'
;	This file can have multiple extensions depending
;	on how many histograms SLUG output. The 0th extension
;	like all structures is empty. The first through nth
;	extensions are the CMDs. Each structure has the following
;	format:
;   TIME            DOUBLE    Array[NSTARS]
;   X               DOUBLE    Array[NSTARS]
;   Y               DOUBLE    Array[NSTARS]
;   X_AXIS_NAME     STRING    'filter1-filter2'
;   Y_AXIS_NAME     STRING    'filter2'
;   
;   Units of TIME are years, and units of X and Y are magnitudes
;	(unless bolometric magnitude or Q(H^0) are set then the
;	units are something odd)
;
;   To find all entries for a particular time (for instance the
;	last timestep), do the following:
;	a=mrdfitS(path+'cmd_'+modname+'.fits', 1)
;	sub=a[where(a.time EQ max(a.time))]
;NOTES
;	as is the case for the SLUG CMD histograms,
;	stars less than 1 Msun are not kept track of
;
;	file size can become enormous!
;
;Written by R. da Silva, UCSC, 11-18-11
;-
pro hess2stars, modname, path=path, abslumcut=abslumcut
pc= 3.0860000d18 ;define a pc fo later use
;------
; Deal with defaults
if ~keyword_set(path) then path='$SLUG_DIR/output/'
if ~keyword_set(modname) then begin
	print, 'model name not set'
	return
endif
if ~keyword_set(abslumcut) then abslumcut=3
;------
; open up the SLUG histo file
h0=mrdfitS(path+'histo_'+modname+'.fits', 1)
nhisto=n_elementS(h0)
;loop over each histogram
for hi=0, nhisto-1 do begin
h=h0[hi]
nt=n_elementS(h.time_axis)
nx=n_elements(h.log_x_axis)
ny=n_elementS(h.log_y_axis)
h.histo=reverse(h.histo)
;convery from logs to magnitudes
x=reverse(h.log_x_axis)
x*=-2.5
y0=h.log_y_axis
y=-2.5*(y0-alog10(4*!pi*(10*pc)^2))-48.6
x0=x
y=y##replicate(1, nx)
x=x#replicate(1, ny)
y3=rebin(y, nx, ny, nt)
wh=where(y3 LT abslumcut AND h.histo,ct)
nstars=total(h.histo[wh])
;define prototype structure
str={time:dblarr(nstars), x:dblarr(nstars), y:dblarr(nstars),$
		x_axis_name:h.x_axis_name, y_axis_name:h.y_axis_name}
k=0
;loop over each timestep and populate array
for i=0, n_elements(h.time_axis)-1 do begin
hist=h.histo[*,*,i]
wh=where(hist AND y LT abslumcut, ct)
 for j=0,ct-1 do begin
   n=hist[wh[j]]
   str.time[k:k+(n-1)]=h.time_axis[i]
   str.x[k:k+(n-1)]=x[wh[j]]
   str.y[k:k+(n-1)]=y[wh[j]]
   k+=n
 endfor
endfor
;write output
mwrfits, str, path+'cmd_'+modname+'.fits', create=hi EQ 0
endfor
end
